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Corrections to the zero-temperature Thomas-Fermi description of a dilute interacting condensed 
Bose-Einstein gas confined in an isotropic harmonic trap arise due to the presence of a boundary 
layer near the condensate surface. Within the Bogoliubov approximation, the various contributions 
to the ground-state condensate energy all have terms of order 7? _4 lnii and R~ 4 , w here R is the 
number-dependent dimensionless condensate radius in units of the oscillator length y/%/mujQ. The 
zero-order hydrodynamic density-fluctuation amplitudes are extended beyond the Thomas-Fermi 
• radius through the boundary layer to provide a uniform description throughout all space. The 

first-order correction to the excitation frequencies is shown to be of order R . 
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I. INTRODUCTION 



The recent experimental realization of Bose-Einstein condensation (BEC) of alkali-metal gases in magnetic 
traps has generated great interest in the physics of a confined, interacting, dilute Bose gas. In the Bogoli- 

ubov approximation [Q, which is generally valid at temperatures well below the BEC transition temperature, the 
macroscopic occupation of the ground state far exceeds that of excited states; the condensate is then described by 
q ■ the relevant Gross-Pitaevskii (GP) equation ^|J|]. This equation has been solved numerically for bosons in both 
isotropic and anisotropic |t0| traps, and analytically in the limits of both a nearly ideal gas (weak interparticle 
interactions) and a dilute nonideal gas (strong interparticle interactions) pd| , [T^] . 

Given the resulting condensate wave function 'J, the linearized small normal-mode amplitudes u and v satisfy 
(*f-} i the coupled Bogoliubov equations ||, which contain the condensate density 2 through the potential energy of 
' interaction. An alternative hydrodynamic approach makes use of the fluctuations in the density and the velocity. 
7— I , It has been shown to be wholly equivalent to the Bogoliubov description p^lo||, and it accurately describes the 
low- lying excited states []l6|-^8| , as recent experiments have verified [|l9|,^0) in considerable detail. 

For large particle number, the mean kinetic energy of the condensate is much smaller than both the interaction 
(Hartree) and trap confinement energies. Neglecting the kinetic energy entirely, corresponding to the Thomas-Fermi 
(TF) approximation, provides an accurate description of the condensate in the interior of the cloud. Near the surface of 
the trapped gas, however, the kinetic and external potential energies become comparable, and the TF approximation 
breaks down. Using a boundary-layer theory, Dalfovo et al. |2l| have calculated the kinetic energy as a function of 
the condensate radius; the formally divergent TF kinetic energy is cut off by a boundary layer of thickness S oc i? -4 / 3 , 
'"^3 ' where R is the (large) dimensionless TF condensate radius. 

The present work extends the boundary-layer formalism of Ref. |2l]] to determine the leading-order corrections to the 
TF description of the condensate wave function, the condensate energy, and the low-lying collective modes. Section 
II summarizes the basic formalism and obtains the first correction to the condensate wave function, both in the bulk 
(of order R~ A ) and in the boundary layer (of order 8). In Sec. Ill, we show that this perturbative expansion gives 
rise to correction terms of relative order R~ A In R and i?~ 4 in the normalization of the condensate wave function and 
in the trap and interaction energies of the condensate. A combination with the previously evaluated kinetic energy 
pl| provides the leading correction to the TF total condensate energy. In Sec. IV, we show that the structure of the 
condensate boundary layer plays an essential role in constructing the corresponding boundary layer for the analytical 
hydrodynamic normal modes u&; the resulting hydrodynamic amplitudes for the density and current fluctuations 
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vanish exponentially as r — > oo, as expected from their equivalence to the eigenfunctions of the coupled Bogoliubov 
equations. 

II. CONDENSATE WAVE FUNCTION 

For a dilute interacting inhomogeneous Bose gas in an isotropic trap potential V ex t at zero temperature, the total 
occupation of the excited states is small, and one can apply the Bogoliubov approximation p|j22|. The spatially 
varying condensate wave function ^(r) is then isotropic and satisfies the (stationary) Gross-Pitacvskii || equation: 

(T + V ext + V B -^(r) =0, (1) 

where T = -h 2 V 2 r /2 m is the kinetic energy, and the trap potential V e xt(i") — muj 2 r 2 /2 is taken to be isotropic for 
simplicity. The Hartree energy is the mean energy of interaction of a particle at r with all the other particles, defined 
as Vn(r) — j d?r' U(r — y')\^(y')\ 2 w ^^(r)! 2 , where the last form reflects the short-range two-body interaction 
U(r) « g8( 3 \r). To leading order, the coupling constant g = 47: ah 2 jra is written in terms of the (low-energy) s-wave 
scattering length a to make contact with experiment; in the present work we only consider a > 0, corresponding 
to interparticle repulsion. The chemical potential fi fixes the total number of condensed atoms Nq — J d 3 r\^(r)\ 2 
through 

fxN = {T} + {V e ^) + {V n ), (2) 

where the noncondensate contribution to the chemical potential is neglected and (■••) = J d 3 r $f* ■ ■ ■ ^ denotes an 
expectation value in the condensate ground state. 

It is convenient to use the the oscillator length do = Wh/mtJo and the oscillator energy hu>o as the basic dimensional 
units. Thus we introduce the dimensionless length z = r/dg, and the dimensionless parameter 

Vo = -j— (3) 
«o 

that characterizes the strength of the interparticle interactions. Experimental conditions typically give t/o 3> 1 , 
corresponding to a strongly interacting system (but the Bogoliubov approximation still requires a dilute system 
with a large interparticle spacing n -1 / 3 relative to a). In this TF limit, the repulsive interactions expand the condensate 

cloud beyond d to a large dimensionless radius R cx A^ J (expressed in units of d ) The kinetic energy in 

Eq. (pi) then becomes small, and the condensate density has a simple parabolic form no = l^l 2 ~ g 1 ^ — V cx t) — 
(/i/<?)(l — z 2 / R 2 ). Since -qo cx Nq scales like R 5 , it is convenient to define 

r/o = r/o i? 5 , (4) 

where fjo approaches a constant value for large R. Defining a scaled variable x = z/R, the full Gross-Pitaevskii 
equation (||) becomes: 

[_ ! e V^ + i* 2 + |%)| 2 -£]#(z) = 0, (5) 

where e = 1/R 4 is a small coefficient in the present limit R — > oo and fl = /.i/R 2 is the scaled chemical potential. 
Here, the scaled condensate wave function 

|g,2 _ ^dfap 2 = 4n4#% 2 

1 1 N Q R 2 1 1 N Q 1 1 W 

becomes independent of R for large R. 

The normalization of the condensate wave function is then written 

/>oo 

r/ ^-^-y o dxx 2 \*(x)\ 2 , (7) 

which defines the condensate radius in terms of the particle number. Correspondingly, the scaled chemical potential 
follows ifrom Eq. @) as p, = (V cxt + Vh +T), where 

(y cxt )^^ = ^J™ dxx^(x)\ 2 - (8) 



2 



_ (Vk) _ N 



R 2 



= — / dx 

Vo 



tf(a;) 



(9) 



its = SI - ^2! 

-R 2 2% 



dxx 



(10) 



where all the energies are expressed in units of Tiluq, and = /dx is the scaled radial derivative. 

In the Thomas-Fermi (TF) limit, one sets the small parameter e to zero in Eq. (|b|). The approximate condensate 

wave function ^>(x) w <&tf(:e)0(1 — a?) = y|(^~ £ 2 )6>(1 — a;) accurately describes the condensate in the bulk Jl0| , pl| , 

where we take p, — ^, defining the condensate radius by = y/2fi. With this choice, the normalization condition 
Eq. (Q) then determines the condensate number N Q in terms of the radius R and the chemical potential. 

The Thomas-Fermi approximation fails near the surface at x = 1, however, where V ex t = \x 2 is comparable to jl; 
in this region, the kinetic term in Eq. (|5|) becomes significant, giving rise to a logarithmic divergence in Eq. (|l0|). For 
small positive e, a boundary layer of thickness 8 forms in the vicinity of x = 1 where the condensate wave function 
varies rapidly. Using standard techniques in boundary-layer theory p3[ , we define an outer solution Pouter = x( x ) 
valid in the bulk region < x < xq < 1, and an inner solution dinner, valid throughout the surface region xq < x < oo; 
an asymptotic analysis matches these two solutions near the boundary x ~ xq. 

The outer (bulk) solution \ may be expressed as a perturbation series in powers of e: 



X(x) = Xo(x) +ex\{x) + ■■ 
Substituting Eq. @ into Eq. (§) yields 

O(e ): X l{x) = \{l-x 2 ); 



< x < xq- 



Oie 1 ): 2 X o(x)xi(x) = 



2 Xo (x 2 - l) 2 ' 



(11) 

(12) 
(13) 



where we have set jl = ^ and specialized to the present case of a real isotropic condensate wave function. In order to 
determine the asymptotic behavior of the outer solution near the boundary, one may write x = 1 + SX with S\X\ <C 1 
and \X\ 1 (where X is large and negative). As X — > — oo, a straightforward calculation reveals 



X (x) ~ *V2 ( _ x) i/2 (i + JL) 



£3/2(_ X )3/2 



21 

8A^ 



(14) 



where the leading-order behavior (oc 8 1 / 2 ) agrees with previous calculations E4). 

The asymptotic behavior of the outer solution implies that the inner solution has the form ^\ ancr {X) = (5 1 / 2 $(X); 
it may be expanded as a series in <5: 



C {X)=8 1 ' 2 [$o(X) + 5$ 1 (X) + 



(15) 



where Xq < X < oo and Xq is large and negative. Substituting Eq. ( [To] ) into Eq. (|^) immediately gives the 
"distinguished limit" |^3| e = 25 3 that balances the leading-order (gradient and nonlinear) terms in the resultant 
differential equation. We readily obtain: 



0(6°) : $S(I)-[I + ^(I)]$ (I)=0; 

0(£ x ) : $'/(A) - [X + 3*g(JT)] $i(X) = -2$£,(X) + ±A 2 $ . 

It is straightforward to verify that, for X — ► — oo, the inner functions ^(A) exactly reproduce Eq. (fl4|): 



(-x) 3 / 2 



i 



21 



8A^ 



(16) 
(17) 



(18) 
(19) 



For large positive X, ^i nneT (X) is asymptotically proportional to an Airy function and therefore decays exponentially. 
Indeed, since Eq. (|l^) defines what is known as the second Painleve transcendent, the boundary condition <I>o(A — > 
oo) V2Ai(A) ~ (27r) -1 / 2 A -1 / 4 exp(— |A 3 / 2 ) must be imposed in order to ensure an unbounded solution for $o(A) 
with no critical points over all A p5[. 
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III. CONDENSATE ENERGIES 



The normalization of the condensate wave function can now be determined explicitly from Eq. (Q) by separating 
the integral into two parts at Xq . An expansion of each contribution shows that the resulting sum is independent of 
the matching point xq = 1 + SXq and yields 



Vo 



1 

15 



5 2 I 



1 



-e In e + e 



1 

24 



In 2 
~~2~ 



J 



(20) 



where / and J are definite integrals involving the inner function &q(X) (see Appendix [A| for details). In fact, / can 
be shown to vanish identically, so that the leading corrections are of order R~ 4 lni? and R~ 4 (instead of order i? -8 / 3 
implicit in the term of order S 2 ). 

A similar analysis yields the physical quantities in Eqs. (p|)-(^0|): 



(Kxt) « §- 



1 „, 1 / 5 In 2 

35 +5/+ 24 £lne + £ Ui" — 



K 



(21) 



(Vh) 



No 

m 



2 
105 



12 



-e In e + e 



In 2 
~3~ 



(22) 



24 £ln£ + e 



1 

12 



In 2 
~6~ 



M 

T 



(23) 



where the constants K, L, and M are also definite integrals involving the inner function <i>o(A). A combination of 
analytical and numerical techniques gives the explicit results 

(24) 
(25) 
(26) 
(27) 
(28) 



/ = 


0; 


J = 


3 t 5 . 

4^ 24' 


K = 


1 T 7 • 

4^ 24' 


L w 


0.4539; 


M = 


1 T 1 

2 12 ' 



such that K + L = J + M. 

Inserting the explicit expressions from Eqs. j2l])-( p3[) in to Eq. (| 
the internal consistency of the calculation. Equation (|20|) implies 



vo(R) = [ Tq )n 



15 



R 
~2 



we immediately recover [i = \R 2 , demonstrating 



^-|ln(1.4128E), 



(29) 



where the constant A is 



A = ±exp(§L 



1.268. 



(30) 



This equation (|29|) relates the condensate number to the chemical potential (and hence the condensate radius); 
conversely, its inverse 



R(Vo) w (1577o) 1/5 + Ml5 Vo r 3/5 In (84.46%) 



(31) 



relates the radius of the cloudfand the chemical potential) to the condensate number. In each case, the first terms 
correspond to the TF result plf . 

The various contributions to the total energy [Eqs. (p|)-(|l0|)] yield 



iV 



3i? 2 
2R 2 



5 , (R 

1 H r m -t 

3R 4 \ A 



R 4 \ A 



(32) 



(33) 



4 



(T) 5 , (R 



AT '~2^ ln UJ' (34) 

and the energies are expressed in units of Hujq. 

The expression for the kinetic energy in Eq. ( |3~I| ) is a small correction that involves only the leading contributions to 
the bulk condensate wave function and the boundary layer; it reproduces the result of Dalfovo et al. pl| ] , and is similar 
to that found in |24| . It should be kept in mind, however, that the present calculation defines the condensate radius 
R in terms of the chemical potential /i = \R 2 \ it therefore includes correction terms beyond the TF approximation, 
as shown in Eq. (|3l|). 

The remaining contributions Eqs. (|3^) and ( J33| ) explicitly require the leading corrections xi and 3>i and have not 
been evaluated previously. The dimensionless total energy per particle {E)/Nq = (T + V cx t + hVn) is found to be 



(E) 5i? 2 I" 4 f R . 



where the first and second terms are, respectively, the standard TF result and the combined contribution from both 
the boundary layer and the first correction to the bulk condensate wave function. Equations ( p9|) and ([35]) together 
constitute a parametric representation of the relation (EJNq)) |2q], a nd it is not difficult to verify that /i = d(E)/dN . 

In fact, Stringari HQ has noted that Eqs. @-(||) and(§|) can be obtained solely from the expression for 
the kinetic energy (|34D . The virial theorem implies that the three contributions to the total energy must satisfy the 
condition 

2<T) + !<Vh>- 2 <^ cxt ) =0. (36) 

The leading-order contributions to the potential energies (14 x t) and (Vh) are known from TF theory, and the correction 
terms must have the same i?-dependence as the kinetic energy but with unknown coefficients. Equation ( |36| ) and the 
condition fi — d{E) / BNq together determine for the unknown coefficients, reproducing the expressions found above 
by direct integration. 

As shown in Fig. |l|, the results of the boundary-layer theory agree strikingly with exact results obtained by in- 
tegrating the GP equation (0) numerically. The analytical slow Rr 2 \nR decay of the average kinetic energy per 
particle (shown dashed bold) with the universal scaling parameter 770 = N^a/do accurately captures the behavior 
obtained numerically (shown solid bold). As a result, the boundary-layer theory provides a much better estimate of 
the total energy per particle (and therefore the chemical potential) than does the TF approximation. In spite of the 
slow decrease in the magnitude of (T)/N with the number of particles, however, it should be emphasized that the 
TF approximation to the total energy and the chemical potential is correct to better than 1% when 770 ~ 1000, due 
to their R 2 increase. 



IV. EXCITED-STATE WAVE FUNCTIONS 



In the Bogoliubov approximation at zero temperature, the interparticle repulsions excite only a small fraction of 
all the particles out of the ground state into self-consistent normal modes. The resulting eigenfunctions u ? - (x) and 
Vj(x) and eigenvalues Ej for the noncondensate modes satisfy the coupled linear Bogoliubov equations HM,E2]. For 
the present purpose, it is convenient to factor out the exact (real) GP condensate wave function 'J, defining sum and 
difference amplitudes |i~5|] 



F^ ^^ and Gj(x)s ^(*)-yW , (37 ) 
*(x) J *(x) 



which are essentially the hydrodynamic amplitudes. Specifically, the perturbation in the velocity potential <pj is 
proportional to Fj (so that the perturbation in the current density is proportional to ^VFj) and the density 
perturbation pj is proportional to ^ 2 Gj. In the TF limit, it is convenient to rescale these amplitudes, with Fj = RFj 
and Gj = Gj / R, yielding the coupled linear equations 

- iV-^VFj) = E^ 2 Gj, (38) 
2^> 4 Gj - ieV-^VGj) = Ej^ 2 Fj, (39) 

where ^ 2 is the (real) scaled condensate density from Eq. (||). These Eqs. (^||) and ([39|) are self-adjoint with a 
normalization integral / d 3 x ^> 2 (F*Gj +G*Fj) — 1 chosen to ensure positive energies Ej for the stable solutions p2[ . 
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In the TF regime (R — * oo), Stringari JT6| has solved these coupled equations for a stationary isotropic condensate, 
but the resulting hydrodynamic eigenfunctions are defined only inside the condensate (x < 1). In contrast, the original 
Bogoliubov equations are well-defined throughout the whole region, including the classically forbidden region (x ^> 1). 
Thus it is interesting to develop a boundary-layer description of the excited states similar to that for the condensate. 



A. Outer (bulk) region 

In the "outer" region (0 < x < xq < 1), we follow the procedure for the condensate and expand the outer scaled 
amplitudes (B7j) and eigenvalues in powers of e: 



■eFl 



Gj ~ G^ ~h cG a 



E 4 



EV + eE* 



(40) 

(41) 
(42) 



Substituting these expressions into Eqs. 
Eq. (|lT|), one obtains: 



) and (B9J) , and including the outer expansion of the condensate given in 



O(e ) 



-V- 



(xlVF°)-(E°) 2 F° = 



j 2 X l ' 



-v-(x^)-(^) 2 ^ 



F° 



1 

"2^ 



2 

Xo 



+ 2V- (xoXi VF° 
£°Xi^° 



2£°^F?; 



Xo 



(43) 

(44) 
(45) 



The spherical symmetry of the condensate density permits a decomposition into spherical harmonics, with the 
normal- mode amplitudes proportional to Yi m (6,<fr). Equations (f43| ) may then be solved explicitly, and the corre- 
sponding (unnormalized) zeroth-order radial amplitudes can be taken as 



F^(x) = p° nl (x)=x l P nl (x 2 ) 



G 



I it 



E nl F nl 

2X§ ' 



(46) 



where the radial quantum number n denotes the number of nodes. Here the eigenvalues take the well-known result 



E 



y/l + n(2n + 2l + 3), 



(47) 



and P n i(x 2 ) — F{~n 1 n + I + |; I + f ; x 2 ) is a hypergeometric function [£8| that terminates to give a polynomial of 
order x 2n . 

While analytical solutions to the inhomogeneous differential equations (^4|) and pq ) are not easy to find, the 
asymptotic behavior of all the outer solutions in the vicinity of the boundary layer x <~ 1 may be readily ascertained. 
Setting x = 1 + SX with X <C — 1 and expanding through 0(S 3 ), we obtain: 



O(e 



0(6*) 



F 



III 



G 



III 



E nl 

25X 


'G° a + 6XG° b + 


F e 
a -+ 


F b F c 

+ + 


S 3 X 3 


S 2 X 2 SX 




f G e a G% 


2SX 


\S 3 X 3 S 2 X 2 



^ + F*ln(-5X) + F^(-SX); 
G 



-±+G% + G% H-5X) + G) \n\-SX) - 



25X ' 



(48) 
(49) 

(50) 
(51) 



where the constants F*_f and are relegated to Appendix [b] for clarity. The expansion (48) is the conventional 

Taylor series of the zeroth-order amplitudes F^(x) about x = 1, where a prime denotes a derivative evaluated at 
x = 1. A straightforward calculation following from the properties of the hypergeometric function 15 2l| shows that 



G 



= (-!)" + 5 A'(D = M J A(1) .... (52) 

with higher-order values listed in Appendix^. Since e = 2(5 3 , the outer solutions in the region x ~ 1 are respectively 
F nl » F° + 2<5 3 i^ and G ni « G°, + 2<5 3 G^ . 

B. Inner (boundary-layer) region 

The need for a boundary-layer description of the eigenfunctions is clear from the form of condensate density 
^(x) 2 — (5$(A) 2 and the gradient Va; = i _1 Vx (arising from the substitution x = 1 + SX). As a result, the second 
term on the left-hand side of Eq. ( |39| ) is now of order e x S/S 2 ~ <5 2 , which is wholly comparable with the first term 
and thus no longer negligible in zero order. It is not difficult to see that the proper scaling of the two inner amplitudes 
in the boundary layer is 

F nl (x) = A nl (X) and G nl (x) = MjQ. (53) 

The gradients in Eqs. ( |38| ) and ( |39| ) must be expanded in ascending powers of 5, leading to the following coupled 
equations for the inner solutions 

d ( ndA n A „ 2$ 2 dA nl ,,/(/ + 1)$ 2 , 



dAV dA y 1 + £A dA (1 + <5A) 5 

d ( dB n ,\ „ 2$ 2 dB nl A(l + 1)^ 2 , 9< 



V y 1 + OT dX (1 + OT) 2 

The asymptotic expressions for the outer amplitudes in the boundary region, Eqs. ([l8])-([5l]), indicate that a simple 
expansion of the inner functions in powers of 6 is insufficient. In order to ensure a match in the vicinity of the 
boundary layer to order 0(e) = 0(S 3 ), the inner solutions must also include the non-trivial terms S 3 \nS and S 3 In 2 5. 
The eigenfunctions and eigenvalues are therefore expanded in the form: 

A nl (X) = A° nl (X) + SAUX) + S 2 A 2 d (X) + 5 3 A 3 nl {X) + S 3 lnM^(A) + S 3 In 2 SA^X), (56) 
B nl (X) = B° nl (X) + 5BUX) + 8 2 B 2 nl (X) + S'B^X) + S 3 hi6B^{X) + S 3 In 2 6B 5 nl (X), (57) 
E nl = E° nl + 5E l nl + 5 2 E 2 nl + S 3 E 3 nl + 5 3 In 8E A nl + S 3 In 2 SE 5 nl . (58) 



Inserting these expansions into Eqs. ( p4| ) and (|5^) , and taking into account the inner expansion for the condensate 
$(A) = $o(A) + 5$>i(X) + ■ ■ •, one obtains equations for the lowest-order inner functions: 

0(6°) : -($ 2 A»/)' = 0; (59) 

- + 2*X* = Kl^Ku (60) 

0(6') : - ($ 2 ^/)' - 2 (® ®iA° nl ')' - 2$ 2 A°/ = 2JS°,$gB° i; (61) 

- ~ 2 ($0*1^/)' - 2<f 2 i?°/ + 2*gfl£, + ml^B° nl 

= E° nl (®lA l nl + 2$o*iJ&) + (62) 

where a prime denotes a derivative with respect to A. 

Equation (|5^) can be solved by taking A^[ as a constant, and comparison with the first term of Eq. ( ff8| ) shows that 
^ni = Pni (I)- ^ n contrast, Eq. (^0|) is inhomogeneous, with the explicit solution 

B o m _ KMV d W) t o o dlnj> (A) 2 



By inserting the asymptotic expression (JTSJ) for <E>o(A), it is simple to verify that B^(X) matches the leading term 
ofEq. @. 



7 



The solution to Eq. ( |6l| ) is readily found to be A^(X) = X, which is identical to the second correction term 

of the outer solution. The other correction B^: d (X) satisfies a somewhat more complicated inhomogeneous equation 
(|62|). While a closed- form solution cannot be found, the behavior in the overlap region X — > — oo is easily obtained: 

BUX) ~ (V - i) - (64) 

Matching the inner solution £?* z with the outer function to order S requires Eh = 0. It then becomes possible to 
write an expression for Bh(X) valid for all X: 

B nl (X)- ^ , (65) 

where the two functions S and Hj satisfy inhomogeneous equations similar to Eq. ( |Tt[ ) for $i: 

- - + (X + 3* 3 )E„ . 4*.*;*, - 2*. A ( X + |l) -A , (66) 



-3;' + (X + 34j)Si = Xt . (67) 

For large negative X, the solutions approach S a ~ i(— X) 1 / 2 and Sb ~ — i(— X) 1 / 2 , so that indeed matches the 
first correction term in Eq. (|49|). Each additional term in 5 may be analyzed in a similar fashion; the procedure is 
straightforward but extremely tedious, so explicit solutions for the inner functions in the overlap region (which are 
lengthy) are omitted for brevity. It is important to note, however, that both inner functions A n i and B n i be properly 
matched to their outer region counterparts at each stage, in order to yield conditions on both the unknown coefficients 
and the eigenvalue corrections. 

By following the above prescription in turn for each term in the inner expansion, it can be shown that all corrections 
to the eigenvalue with prefactors smaller than 0(S 3 ) (including the 0(8 3 In S) and 0(S 3 In 2 5) terms) must vanish 
identically. This result may be formally understood as follows. The density perturbation and the velocity-potential 
perturbation each have corrections to all orders both in the inner and outer regions; in contrast, the only correction 
term in the outer perturbation expansion of the energy ( pl^ ) is of order <5 3 . In order to ensure a smooth asymptotic 
match between the bulk and surface amplitudes to each successive order, all the energy corrections appearing in the 
inner expansion (^8|) with coefficients smaller than S 3 must match to zero (the order S case was considered explicitly 
above). The introduction of logarithmic terms to the outer perturbation expansions would give rise to additional 
contributions in the overlap region, leading to inconsistencies in the asymptotic match. The eigenvalue correction of 
order S 3 is finite, however: E^ — 2E e nl . Thus, 

Em = E° nl + |f (68) 

Eqs. ( |3l| ) and ( |68| ) together yield the number-dependence of the excitation frequencies in the TF limit. The asymptotic 
match reveals an energy correction of order e = i? -4 exists, but unfortunately it does not indicate its magnitude nor 
its variation with n, I. 

The present boundary-layer theory indicates that the energy correction of order 6 3 In S vanishes identically. In 
contrast, a sum-rule approach |T^J23| does yield a logarithmic correction to the eigenvalues, proportional to the 
ratio of the average kinetic energy ( |34| ) and external potential energy (^). This latter approach assumes that a single 
frequency exhausts the sum-rule. Such an assumption is thought to be valid in the hydrodynamic regime where a given 
perturbation excites essentially all of the atoms into a particular low-energy collective mode [^0|,glj. In the vicinity 
of the boundary layer, however, the density of atoms decreases considerably, and the single-mode approximation may 
become insufficient. In practice, any logarithmic correction to the energy eigenvalues would be experimentally or 
numerically detectable only if the magnitude of E^ were strongly dependent upon n and I, or if the number of atoms 
were very low (small R). One may estimate the difference between n = energies obtained with and without a 
logarithmic correction, by assuming E^ » 1(1 — l)/?;/2 jl6| where f3i — logi? and /3/ = 1, respectively. The deviation 
between the two approximations is independent of angular momentum, is largest when 770 ~ 30, and is at most 2% of 
the mode frequency when I = 10. At present, therefore, the data are consistent with either theory. 

The (n, ^-dependence of the energy correction E^ in Eq. ( |S8[ ) can not be found by conventional perturbation 
theory. This situation arises because the explicit integrals, which eliminate the logarithmic divergences, contain the 
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nominally unperturbed inner functions; these functions, in turn, are themselves solutions of differential equations that 
implicitly include the perturbing terms. Consequently, we have considered the readily derived variational expression 



JdV [i* 2 VF; • VFj + 2^G*G 3 + ie* 2 VG* • VG 3 ] 



/ dV "if 2 (G* Fj 



(69) 



which is stationary for small variations about the exact solutions Fj and Gj. As a trial function, we use the unperturbed 
outer solutions in Eqs. (|46]), and the corresponding two-term inner solutions given in Eqs. (56) and ([57]); taken together, 
these expressions constitute uniform unperturbed solutions throughout all space. If the various integrals are divided 
at a point xq = 1 + (5Xo,with Xq -C — 1 and 6\Xq\ <C 1, the correction terms arising from the behavior of the 
outer solution near the TF boundary precisely cancel with those from the inner boundary-layer solutions, leaving 
dimensionless integrals of order In 5 through S 3 In 2 S. The asymptotic match requires that only the term proportional 
to S 3 remains finite, and a detailed evaluation shows that the integrals up to order S 2 InS indeed vanish. The explicit 
calculation to order <5 3 is prohibitive, however, due to the profusion of relevant terms. 

The boundary-layer solutions given above now suffice to determine the approximate density fluctuation amplitude 
p n i oc ^ 2 G n i to order S throughout the whole physical region. Inside the condensate, away from the boundary 
(0 < x < xq < 1), the outer solution is simply the zero-order polynomial p n i( x ) found by Stringari []l6||. The 



uniformly matching inner solution <&{X) 2 B n i(X) in the interval Xo < X < oo (where xq 
and X <C —1) reduces to 



1 + 5X , with S\X Q \ < 1 



Pnl 



P°m(l) + 25p nl '(l)<S> Q (X)Z b (X) + 2Sp° nl (l) 



dX 



<l> (X)Z a (X) - 2$ 1 (X) 



dX 



(70) 



It is not difficult to verify that the third term in Eq. ((70) vanishes for X — > ±oo, although it is nonzero inside the 
boundary region. For large positive X, each term in Eq. (7^) vanishes exponentially, so that the density fluctuations 
become negligible beyond the surface region. 

Figure ^| shows the spatial variation of the hydrodynamic amplitudes in the inner (boundary-layer) region, including 
corrections of order S and S 2 . The velocity potential perturbation (j) n i{X) oc A n i(X) (dot-dashed line) to order 5 or 
S 2 is a linear or quadratic function of the inner coordinate X and therefore diverges at large positive X. In contrast, 
both the density fluctuation p n i(X) oc <& 2 {X)B n i{X) from Eq. ( f7(j| ) (solid line) and the perturbation in the current 
density j n i oc Q 2 (X)A nl (X) (short dashed line) vanish exponentially in the limit X — > oo. In fact, such behavior 
of the density and current-density amplitudes for X — > oo holds to all orders in S as a direct consequence of the 
condensate wave function's exponential decay. The results of the boundary-layer theory differ from those obtained 
within the TF approximation, where the density- fluctuation amplitudes are merely finite at the TF radius. 

The inner solutions provide a more detailed view of the behavior of these amplitudes near the condensate surface. 
In particular, the velocity potential and density perturbation coincide only in the outer region X <C — 1, reflecting the 
fact that in the limit e — > the zeroth-order amplitudes and X^G nl obey the same differential equation (|||). Since 
these unperturbed outer amplitudes are polynomials of order 2n + I in the variable x = 1 + 5X , the inner functions 
must be expanded to at least order fi 2n+l in order to ensure a perfect asymptotic match. 



V. DISCUSSION 



The zero-temperature Thomas-Fermi description of an interacting dilute Bose-Einstein gas confined in an isotropic 
harmonic trap has been extended to include contributions from the condensate surface. For the condensate wave 
function, we have generalized the boundary-layer formalism of Dalfovo et al. |2l[| , obtaining analytic expressions for 
the expectation values of the trap and interaction energy that include the leading corrections due to the surface layer 
and to the bulk condensate wave function. The resulting total ground state energy, which includes all terms of order 
i?~ 4 lni? and R~ 4 , has not been evaluated previously. 

The Bogoliubov equations for the excited states are rewritten in hydrodynamic form and solved to incorporate 
the boundary layer to third order in the boundary-layer thickness S oc i?~ 4 / 3 . This analysis provides a uniform 
extension of the hydrodynamic normal modes found by Stringari fl6(| beyond the TF condensate throughout the 
classically forbidden region. The lowest-order correction to the excitation frequencies has the form E nl /R (namely, 
of order e = R~ 4 ). Although a detailed calculation of the coefficient E nl appears prohibitive, the shift in excitation 
frequencies due to finite- number effects should be relevant for current experiments even when 770 S> 1; both the 
sum- rule approach jl6| and numerical calculations |?2| indicate that E nl increases dramatically with both n and I. 
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APPENDIX A: CONSTANTS FOR THE CONDENSATE 

The constants I-M appearing in Eqs. (|7|)-([l0|) are defined as follows: 

/0 />oo 
-oo ^0 



J 

K 
L 



dX 2$ 



bX$l X 



2 



— - — + / dX 2m 



— oo 

1 



> 



2 8X 

3X 2 1 

2 ~ 8X 



5X$l 



dX 2$ 



7X$t 



dX\^-X 2 - — ) + I dX<S> 



-0) 



M = J 1 dX U$ + X<$> 2 ~^-)+ J™ dX (*< + 



Note that K + L = J + M. The integral / can be evaluated analytically by setting the lower limit to Xo 
integrating by parts: 

X 2 f°° 
1= lim — 2.-2/ rfll$ <- 

Xo^-oo 2 J Xo 

Multiplying Eq. (|l^) by $ and integrating, one readily obtains: 



1 

dXX* & = - 



Xo 



The asymptotic behaviors ( |18| ) and &o(X — ► oo) = immediately give the result 1 = 0. 
The derivation of the expressions for J and K makes use of the relation 



dX $0*1 



23 3 
16 ~ 4 



dXX 2 



dX 2$ 



(Al) 
(A2) 
(A3) 
(A4) 

(A5) 
-oo and 

(A6) 
(A7) 



(A8) 



which may be verified by integrating by parts, comparing with Eqs. (|l6|)-(|l9|), and employing the readily proved 
identity 



/OO -| />oo 

dX&*=-- dX + 
-oo ^ J —oo 



Furthermore, with the identity 



dXX<S>l 



dXX 2 - 



(A9) 



(A10) 



which may b e ea s ily c onfirmed by integrating by parts and making use of the governing equation ( |16| ) for <&o(X), the 
expressions (|A2|)-( A5) are found to be related: 
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J 

K 

M 



3t 5_. 

4^ 24' 

±r - X- 

4^ 24 ' 

±r - -L 

2 12 ' 



(All) 
(A12) 
(A13) 



APPENDIX B: CONSTANTS FOR THE EXCITATIONS 

The constants appearing in the asymptotic expansions of the outer amplitudes F®f and G°' ; e , Eqs. (flS|)-(|5l]), are 
explicitly written (where po(l) = P°z(l) an d E = E^): 



p^(i) = £; 2 po(i); (Bi) 

Po(l) - | [El - ZEl + 1(1 + 1)] p„(l); (B2) 

= | [#o - 10^o + 5£ 2 (ia + 1) + 5) - 13J(I + 1)] p (l); (B3) 

G° a =Po(l); (B4) 

G° b = [El - §] p (l); (B5) 

G c - \ [4 - 5£ 2 + Kl + 1) + 1] (B6) 

= ^ [2Et - 29JSg + 5£ 2 (21(1 + 1) + 19) - 35Z(Z + 1) - 9] p (l); (B7) 

K = 0; (B8) 

F fc £ = 0; (B9) 

F c = 11 [3d - + 31(1 + 1)] po(l); (BIO) 

n=c«\ (bii) 

^ = c 6 e ; (Bi2) 

F / = -11 [4 - El (41(1 + 1) + 3) + 51(1 + 1)] p (l); (B13) 

^ = -|po(l); (B14) 

G l = "I [^o " 12] Po(l); (B15) 

G ' = ~s [^o - 23£; o - ul (l + 1) + 66] (B16) 

°d = ~m N) ~ 15E o - El (31(1 + 1) - 70) + 81(1 + 1) - 66] p (l) + C e a ; (B17) 

G\ = Ct- (B18) 

G / = -il [Eo ~ #o (4i(i + 1) + 3) + 51(1 + 1)] po(l), (B19) 



where C% and G^ are constants of integration. 



[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995). 

[2] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, 

Phys. Rev. Lett. 75, 3969 (1995). 
[3] C. C. Bradley, C. A. Sackett, and R. G. Hulet, Phys. Rev. Lett. 78, 985 (1997). 
[4] N. N. Bogoliubov, J. Phys. (Moscow) 11, 23 (1947). 
[5] E. P. Gross, Nuovo Cimento 20, 454 (1961). 

[6] L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 (1961) [Sov. Phys. JETP 13, 451 (1961)]. 
[7] M. Edwards and K. Burnett, Phys. Rev. A 51, 1382 (1995). 

[8] P. A. Ruprecht, M. J. Holland, K. Burnett, and M. Edwards, Phys. Rev. A 51, 4704 (1995). 

[9] M. Edwards, P. A. Ruprecht, K. Burnett, and C. W. Clark, Phys. Rev. Lett. 77, 1671 (1996). 
[10] F. Dalfovo and S. Stringari, Phys. Rev. A 53, 2477 (1996). 
[11] G. Baym and C. J. Pethick, Phys. Rev. Lett. 76, 6 (1996). 
[12] A. L. Fetter, J. Low Temp. Phys. 106, 643 (1997). 



11 



A. L. Fetter, Phys. Rev. A 53, 4245 (1996). 
W.-C. Wu and A. Griffin, Phys. Rev. A 54, 4204 (1996). 
A. L. Fetter and D. Rokhsar, Phys. Rev. A 57, 1191 (1998). 
S. Stringari, Phys. Rev. Lett. 77, 2360 (1996). 

M. Edwards, R. J. Dodd, C. W. Clark, and K. Burnett, J. Res. NIST 101, 553 (1996). 
P. A. Ruprecht, M. Edwards, K. Burnett, and C. W. Clark, Phys. Rev. A 54, 4178 (1996). 

D. S. Jin, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 77, 420 (1996). 

M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. M. Kurn, D. S. Durfee, C. G. Townsend, and W. Ketterle, 
Phys. Rev. Lett. 77, 988 (1996). 

F. Dalfovo, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 54, 4213 (1996). 
A. L. Fetter, Ann. Phys. (N.Y.) 70, 67 (1972). 

See, for example, C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers (McGraw- 
Hill, NY, 1978), Chap. 9. 

E. Lundh, C. J. Pethick, and H. Smith, Phys. Rev. A 55, 2126 (1996). 
M. I. Ablowitz and H. Seg ur, Phys. Rev. Lett 38, 1103 (1977). 



A. Gonzalez and A. Perez, sond-mat/9802025 , have used Pade approximants to interpolate between the TF limit and the 
nearly ideal gas, in both three and two dimensions (note that their definition of R differs from ours). 
S. Stringari, private communication. 

See, for example, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, edited by 
M. Abramowitz and I. A. Ste gun (Nat. Bur. Std s., Washington, DC, 1968), Chap. 15. 



F. Zambelli and S. Stringari, |;ond-mat/9805068 
R. D. Puff, Phys. Rev. 65, A406 (1965). 
S. Stringari, Phys. Rev. B 46, 2974 (1992). 
B. 1. Schneider and D. L. Feder, unpublished. 



FIG. 1. The contribution of the kinetic energy in units of hu>o is shown as a function of the universal scaling parameter 
■qo = Noa/do- The kinetic energy per particle is calculated numerically by direct integration of the GP equation (solid bold) 
and analytically by the boundary-layer theory (dashed bold), whose result is given in Eq. (|34|). The dashed thin line is the 
difference between the exact total energy per particle (computed numerically) and the full boundary-layer theory in Eq. (fjil). 
The solid thin line corresponds to the similar difference keeping only the approximate result of the TF theory [the first term of 
Eq. (53)]. The total energies per particle for a restricted range of rjo are shown in the inset for the exact (solid), boundary-layer 
(dashed) , and TF (long dashed) theories. On this restricted scale, the curves for the boundary-layer and exact results coincide. 



FIG. 2. The (unnormalized) hydrodynamic amplitudes and condensate wave function in the boundary layer are shown as a 
function of inner coordinate X. Bold and thin lines correspond to results calculated numerically to order 8 and 5 2 , respectively. 
The universal parameter is 770 = 1000 giving 5 ~ 0.06, and therefore X — —8 corresponds to x « 0.5. With units chosen 
appropriately, the velocity potential to order S (bold dot-dashed line) and the density perturbation [bold solid line to order S, 
from Eq. (|70j)] coincide in the overlap region X <C — 1; while the former diverges as X — > 00, the latter decays exponentially. 
The perturbations in the current density (dashed lines) to order 8 (bold) and 8 2 (thin) are respectively linear and quadratic 
in the overlap region and decay exponentially at large distances. The results are presented for (n, I) = (0,2); the overall sign 
of the inner hydrodynamic amplitudes is odd in n, and the magnitude of the asymptotic slope (for large negative X) increases 
with I. The inset shows the various amplitudes in the bulk region as a function of outer coordinate x. In the outer region, the 
velocity potential and density perturbation coincide (shown as the solid line) ; the former matches smoothly with that from the 
inner region to order S 2 (thin dot-dashed line). A perfect asymptotic match of the outer current density cx x(l — x 2 ) (dashed 
line) to its inner counterpart would require an inner expansion to order 8 3 . 
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